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Abstract 

The existing 2-D o-/x scheme and o-e scheme based on the method of space-time con- 
servation element and solution element, which were constructed for solving the linear 2-D 
unsteady advection-diffusion equation and unsteady advection equation, respectively, are 
tested here. Also, the a-e scheme is modified here to become the v-e scheme for solving 
the nonlinear 2-D inviscid Burgers equation. Numerical solutions of six test problems are 
presented in comparison with their exact solutions or numerical solutions obtained by tra- 
ditional finite-difference or finite-element methods. It is demonstrated that the 2-D < 2 -/ 2 , 
a-e and v-e schemes can be used to obtain numerical results which are more accurate than 
those based on some of the traditional methods but without using any artificial tuning in 
the computation. Similiax to the previous 1-D test problems, the high accuracy and sim- 
plicity features of the space-time conservation element and solution element method have 
been revealed again in the present 2-D test results. 
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Introduction 


The space- time conservation element and solution element (to be abbreviated as CE/SE) 
method is a method newly developed by Chang[l-3] for solving conservation equations de- 
scribing the flow properties. The concept and methodology in this method are significantly 
different from those in the well-established traditional methods such as finite difference, 
finite volume, finite element and spectral methods. First, the flux is conserved in time and 
space because time and space axe unified and treated equally. Second, all the dependent 
variables and their derivatives are considered as individual unknowns to be solved simul- 
taneously at each grid point. And third, the concepts of conservation element (CE) and 
solution element (SE) are introduced to enforce both the local and global flux conservation 
without using any ad hoc techniques such as interpolation or extrapolation. It has been 
proved that this method is superior to many of the traditional methods in accuracy, sim- 
plicity and generality. Also, the stability conditions for the numerical schemes based on 
this method are the same as or less restrictive than those based on the traditional methods. 
The detailed descriptions of this method can be found in [1-3]. 

All of the 1-D numerical schemes constructed on the basis of CE/SE method have been 
tested to show that they are surperior to some traditional methods in both accuracy and 
stability [2,4]. Recently the two-dimensional versions of the a-/j scheme, a-e scheme and 
Euler solver have also been constructed by Chang [3] based on the same set of design 
principles. 

In this work, the aforementioned 2-D schemes and an extended scheme are tested using 
several model problems. The 2-D unsteady advection-diffusion equation is solved by use 
of the a-fj, scheme. The influence of diffusivity on the accuracy of the numerical solution 
is investigated by comparing with the exact solution. The 2-D pure advection equation 
is solved using the a-e scheme with varying e to study its effect on numerical diffusion. 
Then the 2-D pure diffusion equation is solved using the \i scheme, which is a special 
form of the a-fj, scheme by letting a x = a y = 0. The numerical results are compared with 
those obtained by Crank-Nicolson implicit finite-difference and finite-element methods. 
Furthermore, the a-e scheme is modified to become the v-e scheme for inviscid Burgers 
equation. Two test problems are used to demonstrate that the v-e scheme is still more 
effective than several traditional finite- difference methods in solving nonlinear problems in 
a much simpler manner. 
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Numerical Schemes 


1. The a-fj, Scheme for Advection-Diffusion Equation 

The basic a-fx scheme is first reviewed, which is constructed for solving the linear 
advection-diffusion equation 

du du du d 2 u d 2 u. 

ft +a *te + a «^ = ,x{ d7i + df 2) 

where a x and a y axe the constant advection speeds and fi is the viscosity coefficient. In the 
space-time Euclidean space E 3 , the integral form of (1.1) is 

J J v (V -h)dV = 0 (1.2) 

where h = ( a x u — fxdu/dx , a y u — fidu/dy, u ) and F is an arbitrary space-time region in £3. 
By use of Gauss’ divergence theorem, Eq. (1.2) becomes 

<f h • ds = 0 (1.3) 

/s(V) 

where 5(F) is the boundary of region F and ds = dan with da and n, respectively, being 
the area and the outward unit normal of a surface element on 5(F). 

The conservation element (CE) and solution element (SE) are the two basic elements 
used in the construction of numerical schemes. Some representative CE(j, n) and SE(j, n) 
are depicted in Fig. 1. There are two sets of mesh points (j, k, n ) where j and k are spatial 
mesh indices and n is the time mesh index(see Fig. 1(a)). One of which, denoted as ftj, 
is the set of solid circles with j, k = 0, ±1,±2, ••• and n = ±1/2, ±3/2, ±5/2, • ■ -. The 
other denoted as fl 2 , is the set of harrow circles with j, k = 1/3, 1/3 ±1,1/3 ±2, ••• and 
n = 0,±1,±2, •••. Each mesh point (j,k,n) in fix or fl 2 is associated with three CEs, 
which are quadrilateral cylinders in space-time, one SE, which is the union of three vertical 
planes, a horizontal plane, and their immediate neighborhood. The CEs and SEs in fli, 
denoted as CEp^(j, k, n + 1/2),/ = 1,2,3 and SE M(j,k,n + l/2)(see Fig. l(b)(c)), are 
different from those in fl 2 , denoted as CE^(j, fc, n + 1), / = 1, 2, 3 and SE^ 2 i(j, k,n + l)(see 
Fig. l(d)(e)). The geometry of the hexagon ABCDEF, which appears in the CEs(SEs) 
at the mesh points (j, k, n) € both fix and fl 2 , is determined by three positive parameters 
w,b and h(see Fig. 1(f)). 

For any ( j , k, n ) € ft, let 


OEV- 7 x SEW(j,k,n), if (j,k,n)e(l 1 
SE(j,k,n)=< 

5E^(i,fc,n), if (j, k, n) £ ft 2 


(1.4) 
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For any (x ,y,t) € SE(j, k, n), u(x,y,t) and h(x,y,t ) are approximated by 

u*(x,y,t-,j, k,n ) = + (u x )l k (x - Xj , k ) + (%)£*(</ - ?/;,*) + (u t )^(i - t n ) (1.5) 

and 

h*(x,y,t;j,k,n) = [ a x u*(x, y, t; j, k, n) - ydu*(x, y, t; j, k, «)/ dx, 

y, t; j, &, n) - /*3u*(x, y, t ; j, k, n)/ %, tt*(x, t/, i; ;, A:, n) ] (1.6) 

where u^ k , (u x )% k , (%)£* and (u t )] ik are constant within SE (j,k,n). The assumption that 
u — u*(x,y,t]j,k,n ) satisfies (1.1) implies 

( U t)j,k = ~ \ a x(Ux)j jk + 0,y{u y )j k (1-7) 

Substituting (1.7) into (1.5), one has 


u (x, y, t; j, k, n) - u jtk + (u x ) jyk [(x Xj,k) o, x {t t )] + (%)£fc[(2/ — yj,k) — a y (t — t n )] (1.8) 

Thus there are three independent marching variables, i.e., uj k , (u*)” fc , {u y ) r - k associated 
with a mesh point (j, k, n ) € Cl. 

Considering (1.3) in any arbitrary CEp^j, k,n + 1/2) <E Cli, the approximation is 


Ts(C E?\j, k,n+l/ 2 )) 


h* ■ ds = 0, l = 1,2,3 


Similarly, in any arbitrary CEp^(j, k,n + 1) € Cl 2 , the approximation 


Ts(CEf\j,k,n+l)) 


h*-ds = 0, 1 = 1,2,3 


( 1 . 10 ) 


As a preliminary for the evaluation of the integrals (1.9) and (1.10), a new coordi- 
nate system (CiV) shown in Fig. 1(f) is introduced. Some transformed coefficients in the 
(£, y ) coordinate which will be used in the following computations are described briefly as 
following: 




where 


VJ \a Vl 


w — b w + b 


= T tUx 


' A( ( w + b)A( 

T -1 = 2u> 2 wh 

Ay ( w — b) Ay 

' 2w 2 wh 


( 1 . 11 ) 


( 1 . 12 ) 
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with. 


AC = \J[w — b ) 2 + h 2 , At] — yj(w + fe ) 2 + h 2 
and T l is the transpose of T. Let 

A r\ 


+ AC + 

U t = < 


"Uri 


3 At 


3 Ai 


Uc ~2AC a ° Vr> ~2Ar } ar> 


and 


«c = ^(A 0 2 , ?, = = 


9fj,At 


8w 2 h 2 v w 8uW " ’ ST 8w 2 h 2 

where Ar = 2\/& 2 + A 2 - Furthermore, let cr$ + , <r^ _ , • • • be defined by 


(At) 2 


off* 

— 

1-U C - Ur, 

off* 

= 

±(1 -V^- Z^rj)(l + I'c) + in + Ct — Cc 

a 13 

- = 

d=(l — U{ — Ur,)( 1 + u v ) + Cc + £r — Ct? 


= 

1 + ^ 


= 

T(1 + &0( 2 - t'c) - 2£„ 

4l }± 

= 

±(i + + u n ) + Cc - & + 6, 

°31 

= 

1 + ^7) 

off* 

= 

±(1 + ^„)(i + v ( ) + Cc - + C»? 

off* 

= 

dX 1 + u v )(2 - u v ) - 2Cc 

off* 

= 

1 + u c + u v 

„(2)± 

<r 12 

= 

= F (1 + v <, + v v){^ — u <) — £v — €r + £<: 

off* 

= 

T(1 d* V C d" z/ 7))(l — v v ) — ^ — Cr + Ct; 

off* 

= 

1 -I/ c 

off* 

= 

±(1-^X2 + v ( ) + 2C, 

JVi 

a 23 

= 

=F(1 - I'cXl - ~ & + £r ~ 6, 

off* 

= 

1 -Vr, 

off* 

= 

■f(i — v i j)(f — ^c) ~ £c d" £r — C»? 

off* 

= 

±(1 - z / tj )(2 + + 2 Cc 


(1.13) 


(1.14) 

(1.15) 

(1.16) 


(1.17) 

(1.18) 

(1.19) 

( 1 . 20 ) 
( 1 . 21 ) 
( 1 . 22 ) 

(1.23) 

(1.24) 

(1.25) 

(1.26) 

(1.27) 

(1.28) 

(1.29) 

(1.30) 

(1.31) 

(1.32) 

(1.33) 

(1.34) 
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With the above definitions, the final results of (1.9) and (1.10) are: 


4i )+ “ + a}2 )+ u^ + crg )+ u+]J^ 1/2 = 

K H “ + ^ “< + dJNf]* +1A1M . 1/s 

(1.35) 

<4i )+ « + 42 )+ «< + a 23 + ^]^ 7 = 

[o-W'ti + + ' T S i ~“n] ) _ 2/ ,,i , +1/ 3 

(1.36) 

4l )+tt + <? 32 + U( + crS )+ U+] n ^ X/ = 

j JyK 

K’-O + + 4i ) ‘<]” +1/ 3 t . 2/3 

(1.37) 

[4 2)+ « + &i2 +u t + 4§ )+ *4]”* 1 = 

H?'« + '<] ,*-1,3 

(1.38) 

[4 2)+ « + 4 2)+ «* + 4 2)+ ^]JJ 1 = 

H?'“ + 

(1.39) 

[4l )+u + 032 )+w C + <7 33 + wJ']^ 1 = 

H?'“ + + 4 ) '<]j! 1 1 / / 3 a i+2/ 3 

(1.40) 


Let 4 X) , 4^> 4 2 \ and denote the expressions on the right sides of (1.35)-(1.40), 

respectively, and 

Pn = [3(1 + V <)Q + v v) + (3 + 2 + v v )^ + (1 + Vr,){(,n — £r)]/A (1) (1-41) 

/& = [-3(1 + ^)(1 *0 - (3 - 2«, c - *06 + (1 + *,)(*, - e,)]/A (1) (1.42) 

4? = [(2^ + v v )U + (2- - *,)(& - &)]/A« (1.43) 

4V = [3(1 + i / c)(l + v v ) -f (3 + + 2vr,)Z v + (1 + *0(& — £»•)]/ A^ (1-44) 

Ps 2 = [(v< + 2*)tn + (2' - *0(6- - ^c)]/A (1) (1.45) 

4s = [-3(1 + y<)(l -v < - v v ) - (3 - v c - 2v v )( ri + (1 + - 6)]/A (1) (1.46) 

4? = [-3(1 - *0(1 - v v) - (3 - 2^ c - *,)fc + (1 - Vv)(tr - ^)]/A( 2 ) (1.47) 

P22 = [3(1 — V v ){ 1 + V( + t'n) + (3 + 2v<; + V v )^ + (1 — V v )(Zr — £„)]/ A^ (1-48) 

4s = K 2 ^ + u v)£< + ( 2 + v v )(— & + ^)]/A^ (1-49) 

4? = [-3(1 - *0(1 - v v) ~ (3 - ~ 2^, + (1 - »*)(-& + £r)]/A< 2 ) (1.50) 

4? = [(*< + 2 !/„)£,, + (2 + i^)(— £r + &)]/. A^ (1-51) 

4 2 3 } = [3(1 — ^)(1 + + I'tj) + (3 + + 2 U v )(r, + (1 - v<)Ut — €c)]/ A^ (1-52) 

where 

A (1) = 3 [3(1 + i^)(l + z^Xl — t'c - ^7?) + 2(1 + ^)(1 — 

+2(1 + *„)(! - v t - !/,){, + 2(1 + »<)(! + »„)& + (^) (1.53) 
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(1.54) 


A< 2 > = 3 [3(1 - i/<)(l - i/„)(l + i/ f + v v ) + 2(1 - v ( )(l + u ( + p n )t< 

+2(1 - ^(l + V( + l / v )€r, + 2(1 - I/ C )(l - I+)£r + (^2 w h) 


With the assumption of AW ^ 0 and ^ 0, the final form of the a-p scheme are 


n+1/2 _ 
U j,k ~ 

4 ,) +4 i) +4 1, ]/3 

(1.55) 

Wit 1 * = + $ 4" 

(1.56) 

«)»r /2 = /-SM 1 ’ + <®4" + 484" 

(1.57) 

where (j, k,n + 1/2) £ fii, and 



“"t 1 = [ 

4 2) +4 2) +4 2) ]/ 3 

(1.58) 

(t*?® 1 = 4 2) 4 2) + /®4 2) + 4 2 3 } 4 2) 

(1.59) 

K® 1 = p&* i 2) + 4 2) 4 2) + 4! ) 4 2) 

(1.60) 


where ( j , k, n + 1) € ^ 2 - 

The a-p scheme is second-order accurate in space and time. Detailed derivations and 
discussions of the stability condition for the a-p scheme are referred to [3]. 


2. The a-e Scheme for Unsteady Advection Equation 
Equation (1.1) becomes the pure advection equation 


du du 
dt + Qx dx 


du 
' dy 


+ &y~7T~ = 0 


( 2 . 1 ) 


when pi = 0. It is expected that the a-p scheme may become unstable as p — * 0 for the 
nonlinear problems. Some modifications are thus needed for its improvement. The CEs in 
the a-p scheme are replaced by 


C£ (1 >(j, it. n+ 1/2) 

= [CB^y, k, n + 1/2)] U [CE?\j, k,n + 1/2)] U [Cff’y. k, n + 1/2)] (2.2) 

where (j, k,n + 1/2) € Oj, and 

CEC\j,k,n + 1) 

= [CE{ 2 \j, k,n + 1)] U [CEi 2 \j, k,n + 1)] U [CE®{j, k,n + 1)] (2.3) 
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where ( j , k, n + 1) € ^ 2 - The total flux leaving the boundary of any new CE is assumed to 
be zero, i.e., 


i 


S(CEW(j,k,n+ 1/2)) 


h* • ds — 0 


and 


f 

Jsi 


h* ■ ds = 0 


/S(C23(2)0\fc,n+l)) 

Substituting (1.17)-(1.34) into (1.35)-(1.40) with = £„ = £ T = 0, one obtains 


[« + (! + + (1 + ^)<] =[«-(! + - (1 + *„)< 

[“ " (2 - + (1 + ^K]JJ 1/2 = [« + (2 - ^ ( )uf - (1 + v v )u+ 

[« + (! + v c )v% - (2 - 7 =[«-(! + + (2 - ^X] J +1/3>fc _ 2/ 3 

[u - (1 - I/ C )uf - (1 - =[« + (! “ ^cX + (! - ^)<]^ /a 


j+l/3,fc+l/3 


j-2/3,fc+l/3 


in+1 


“ + ( 2 + ^X “ (1 “ ^r7)X] i>A . = [« - (2 + ^cX + (1 “ ^X]" 


j— l/3,fc— 1/3 
ti+1/2 

.7+2/3, fc— 1/3 
I 71+1/2 


(2.4) 

(2.5) 

( 2 . 6 ) 

(2.7) 

( 2 . 8 ) 

(2.9) 

( 2 . 10 ) 

( 2 . 11 ) 


[“ - (1 - , 'tK + + ( 2 + ">/]■( =[« + (!- D”? " (2 + 

when |i/ c | ^ 1,|^| ^ 1 and |i/ c + i/„| ^ 1. Let 4 1) )4 1) .4 1) >4 2) )4 2) > and s'J 2) denote the 
expressions on the right sides of (2.6)— (2.11), respectively. From (2.6)-(2.8) and (2.9)- 
(2.11), respectively, one obtains 


“S ,/2 = | [(i - *< - «,)$' + (i + ■'«)<?> + (l + x,M 

« + )S 1/! = Hf’-# 1 ] 


3 
1 

3 

/,.o+l n +l/2 _ J - 

(D? /i,* - 3 


(1) 

3 


c'(!) 

S 1 


j(l)l 

•S 3 


where (j, k,n + 1/2) € fii, and 

“i.t 1 = 3 [(! + + w n)4 2) + (1 - ^ C )4 2) + (1 - ^)4 2) 

(“DS 1 = i[4 2, -4 2) ] 


(v)S = 


_ I [ 4 ( 2 ) _ 4 ( 2 ) 


( 2)1 

1 


( 2 . 12 ) 

(2.13) 

(2.14) 

(2.15) 

(2.16) 
(2.17) 


where ( j , k, n + 1) € 



For any (j, k, n + 1/2) € fii, the numerical analogues of and (u+)"| 1//2 using 

the central- difference approximation are defined as 


« + )T 2 = (*&&«/. - “S/iW)/ 6 


and 


« + )?,r /2 = (-SJ&w. - «SXv/« 

where "j+t/M+i/3’ “j-?/«+i/3. and “£i/3?lc-2/3 can be com P u,ed usin 8 

'n+ 1/2 ( \ ^ V 

Then the first marching step of the a-e scheme can be expressed by (2.12) and 

(«< + )?,r /2 = « + )s i/2 + (£ - i/ 2 M<i<)”t i/2 

M)”t i/2 = («; + e ,/2 +(e-i/2)(*j)"t i/2 

where 0 < e < 1, which is a parameter for controling the numerical diffusion, and 

t 1/2 = 2 [(f< + )S ,/2 - « + ® 1/2 ] 

(<*ojr /2 = 2 [« + )s ,/2 - K + sr /2 ] 

For any (j,k,n + 1) € fl 2 , the counterparts to (2.18) and (2.19) are 

(“c + )”,r = (“Srm/j - “#/«-./.)/« 

and 

( W tj + )j,fc 1 = ( W j-"l/3,*+2/3 — U ?-l/3,k-l/z)/® 
where u£t/ 3 ,*-i/ 3 > w j+ 2 / 3 ,*-i/ 3 > and M j-t/ 3 ,fc+ 2 / 3 can be obtained by using 

At s n+1/2 


/ Af \ n1 

•tf 1 = (“ + !-<*.)* 


Then the second marching step in the a-e scheme is formed by (2.15) and 


(«?)3‘ = (o?© 1 + (« - 


= (“i + )?r+(£- i/2H<i<)",t 


+ \n+l 


(2.18) 

(2.19) 

(2.20) 


(2.21) 

( 2 . 22 ) 


(2.23) 

(2.24) 

(2.25) 

(2.26) 

(2.27) 


(2.28) 

(2.29) 
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where 


(<*“{■)", t 1 = 2 [ 

K+ffit 1 - («r)3 i ] 

(2.30) 

= 2 [ 

(«; + )?r - w + )",f] 

(2.31) 


Next, the weighted-average counterparts to (u ^ + )” t k and (u'+)j k in the a-t scheme are 
introduced. As a preliminary, for (j, k,n + 1/2) € Cli, let 


Xl 

x 2 

xz 


'n+l/2 

U j+l/3,k+l/3 


— U 


n+l/2 

i,k 


'”+1/2 

2/3,fc+l/3 


— U 


n+l/2 

j,k 


'”+ 1/2 

u j+l/3,k-2/3 


- »/ n+1/2 
U j,k 


and 


where 


»)S 1/2 = + (4°) 2 


I n+l/2 


j,k 


, Z = 1,2,3 


(««)ir /2 = + -a) (»w« 1/2 = 

t..( 2 )\n+l /2 _ 3a?i ( ,( 2 )\n+l /2 _ (36 + U>)®i + 2 u>X 3 

^ - 2w [Uy )j ’ k ~ 2i oh 

s..(3)\n+i/2 _ 3®1 / '..f3hn+i/2 _ (w — 36)®i -)- 2u;a:2 

h,k - 2^ - 2^h 

For any a > 0, the weighted-average counterparts to and u'+ are 


(2.32) 

(2.33) 

(2.34) 

(2.35) 


(2.36) 

(2.37) 

(2.38) 


uJ + = 


°, 

i (fl 2 03H 1)+ + (Ml)"^* + (M2)*t4 3)+ 

(^1^2)" + (^2^3)" + (030l)“ 


if 0 1 = O 2 = ^s = 0 
otherwise 


and 


u™ + = < 


0, 

^z)^ + {036lYuW+ + (^i^ 2 )^( 3 H 
(^1^2)“ + (^2^3)“ + (^3^1)“ 


if 0\ — $2 = @3 — 0 
otherwise 


where 


(“? H )m I/2 = -(2X2 + x 3 )/6 (t.<*> + )”t I/2 = -(X2 + 2 x 3 )/6 

(“f ,+ )?.t 1/2 = (2xi + x 3 )/6 « 2 ) + )”t ,/2 = (xj - x s )/6 


(2.39) 


(2.40) 

(2.41) 

(2.42) 
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(<4 3)+ ffi 1/2 = (** - *>)/6 (« f + )?1 1/S = (2*1 + * j )/6 (2.43) 

For ( j , fc, n + 1) € CI 2 , there are 


_ 7 /»+l _ 7/ n+l 

2/1 — “ 7 — 1 / 3 ,*— 1/3 

(2.44) 

„/n+l „,ra+l 

2/2 — “j+2/3,*— 1/3 “j,* 

(2.45) 

7/o _ 7 >+l _ „»»+l 

2/3 — “j— 1/3, *+2/3 “j,* 

(2.46) 


The counterparts to (2.41)-(2.43) and (2.36)-(2.38) are, respectively, 

(«? )+ )SP = (22/2 + 2/3)/ 6 (t#*® 1 = (2/2 + 2y 3 )/6 

= -(22/1 + y 3 )/6 («i 2)+ )"t 1 = — (yx - y 3 )/6 

(“c 3)+ )”,t 1 = -(2/1 ~ 2/2)/6 («i 3)+ )i,t 1 = — (2yi + ya)/6 


and 


(2.47) 

(2.48) 

(2.49) 


(• W -£<» + ») (2.50) 

(2-51) 

( 4 ?>®' = -£ (W - -^1 a-— < 2 - 52 > 

Equation (2.35) is still valid after changing the index (n + 1/2) into (n + 1). With the above 
definitions, (2.39) and (2.40) can be used as the the weighted-average part for a mesh point 
(j, k, n + 1) € 

Finally, the weighted-average a-e scheme is generalized. The first marching step is 
formed by (2.12), 


(<4)S 1/! = (“DID 2 + (« - i/a)(*4)S v ’ 


and 


«e ,/2 = «)S i/! + ( £ - i/2)(<i«t)it i/2 

where (/, &, n -f 1/2) € Hi. The second marching step is formed by (2.15), 
and 

Wffi 1 = (“Dm 1 + (« - i/2)(<K)Jf 


(2.53) 

(2.54) 

(2.55) 

(2.56) 
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where (j, k,n 4- 1) € 0 2 . 

The accuracy of the a-e scheme is the same as that of the a-fi scheme. The stability 
condition is discussed in [3]. 


3. The fj, Scheme for Unsteady Diffusion Equation 

The unsteady diffusion equation is a typical parabolic second-order partial differential 
equation. It is a special form of Eq. (1.1) with a x = a y = 0, 


du _ . d 2 u d 2 u . 

~dt = + ~8y i ' 


(3.1) 


For the interior CEs, the fj, scheme is a special form of the a-fi scheme by letting 
i/£ = u v = 0. For the boundary CEs, the scheme depends on the boundary conditions. If 
one flow variable is unknown at the boundary mesh points, only one conservation relation 
is needed. Referred to Fig. 2, for the boundary mesh point (j, k,n + 1/2) € ffo, Eqs. (1.35), 
(1.36) and (1.37) can be obtained on the boundary at x = 0 with k = 2,4, 6, • • • ,ny — 1, 
at y — 0 with j = 1, 3, 5, • • • , nx — 1, and at y = 1 with j = 1, 3, 5, • • • , nx — 1, respectively. 
For the boundary mesh point ( j,k,n + 1) € f& 2 , Eqs. (1.38), (1.40) and (1.39) can be 
obtained on the boundary mesh points at x = 1 with k = 3, 5, 7, • • • , ny — 2, at y = 0 with 
j = 2, 4, 6, • • • , nx — 2, and at y = 1 with j = 2, 4, 6, • • • , nx — 2, respectively. 

It is assumed that the value of {u x ) r - k is unknown on the boundaries at x — 0 and x = 1 
and the value of {u y )^ k is unknown on the boundaries at y — 0 and y — 1, respectively. 
And we know 


( u Vlk = [(w - b){u x )l k - h(u y )l k ] /6 (3.2) 

= [(w + b)(u x )l k + h(u y )l k ] /6 (3.3) 


Substituting (3.2) and (3.3) into (1.35)-(1.40), one can obtain 

64 x) - 6o-Jl )+ <t 1/2 + ha{ 2 )+ K)£l 1/2 - haQ + {u y )i+ k 1/2 


( U ®)l ,k ~ 


<7l2 )+ (^' - fy + 43 + (^ + 4 


I S"+l/2 _ 64 1 * ~ 6<rg )+ ii^i /2 - (m - 6 )i4 H (m*)S /2 - (m + t)<rfi )+ (u I )": 

l ’ ,>im m-<4 )+ + 4 ,+ ) 


n+1/2 

ny 


(3.4) 


, ,M-i/2 64 1 ’ - 64 i )+ «jJ 1/2 - ,,, C l 

K)w < 3 - 5 ) 


(3.6) 
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(3.7) 


for the boundary mesh points (j, k,n + 1/2) € Qi, and 


K)£* = 


e4 2) - e <r{? + <s + - MfWS 


, \n+l _ 6 s 3 2) ~ 64 2)+ «i,X 1 - (W - 6)4 2)+ («r)j I l 1 “ (*0 + ^)4s + (^)",1 1 , 0 on 

K)il1 ~ K-*g ) + + 4 )+ ) ( } 

, \n+l 65 2 2) - 6<rg )+ «^S - (™ - %§ )+ («*)£& - (to + 6)^3 + (^)”,tv /« nN 

K)iw ~ m-4 )+ + 4 )+ ) 

for the boundary mesh points ( j , fc, n + 1) € JV 

Equations (3.4)-(3.9) are the // scheme for the boundary mesh points. Once (u x )™ k and 
( u y)j,k become known, and ( u v)],k can be obtained using (3.2) and (3.3). Note that 

z/^ = v n = 0 in oi^ + , cr$ + , • • • and • • • in (3.4)-(3.9). The y scheme is characterized 

only by the diffusion numbers £ r and is unconditionally stable with the same accuracy 
as the a-fj, scheme. 

4. The v-t Scheme for Inviscid Burgers Equation 
The inviscid Burgers equation in conservative form is 

Tti + Tx + Ty =(> (41) 

where f = 0.5u 2 . Its integral form in the space-time Ez is the same as (1.3), except 
h — (/, /, u) here. Some additional definitions cooresponding to / are 

f u = df/du = u, / = 0.5/“ u (4.2) 

fx = df/Bx = f u u x , fy = df/dy = f u Uy , ft = df/dt = /“u t (4.3) 

Similarly, f(x,y,t;j,k,n) is approximated by 

M) = /?* + (/.)7^(* - *;,*) + (fy)Uv - v&) + (/*)£*(* - * n ) ( 4 - 4 ) 


h*(x,y,t;j,k,n ) = (f*{x,y,t;j,k,n),f*(x,y,t]j,k,n),u*(x,y,t]j,k,n)) (4.5) 

Note that (/*)£*, (/ y )**, and (/*)** are functions of (u)**, («*)£*, K)’*, and («*)**. From 
the assumption of 

du*(x,y,t-,j,k,n ) df*{x,y,t;j,k,n) df*(x,y,t]j,k,n) _ ^ 

di da; dy 
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one obtains 


Ut = -fx -fy- ~(f U U x + f U Uy) (4.7) 

One can conclude that there are also three independent variables, namely (u x )^ k , 

and (u y )^ k , associated with each mesh point (j, fc, n) € f l. 

Furthermore, the corresponding coefficients in the ((, 77 ) coordinate are 



Then the counterparts to (1.17)— (1.34) are: 


4* = 1 - 0.5/ c + - 0.5/”+ (4.9) 

4* = ±(l-/ c+ -/” + )(l + / c+ ) (4.10) 

4* = ±(i-f + -r + )(i+f v+ ) (4.n) 

4* = l + 0.5/ c+ (4.12) 

4* = T(l + / C+ )(2-/ c+ ) (4.13) 

4* = ±(1 + / C+ )(l + f v+ ) (4-14) 

4* = 1 + 0.5/”+ (4.15) 

4* = ±(l + /”+)(l + / c +) (4.16) 

4 ± = +(l + /”)(2-/”+) (4.17) 

4 ± = 1 + 0.5/ c + + 0.5/”+ (4.18) 

4* = =F(1 + / C + + r + )(l - / C +) (4.19) 

4* = +(l + / c + + /”+)(l-/”+) (4.20) 

4+ = l-0.5/ c + (4.21) 

4* = ±(l-/ c +)(2 + /+) (4.22) 

4 * = +(l-/ c +)(l-/”+) (4.23) 

4+ = 1-0.5/”+ (4.24) 

4* = T(1 - r + )(l - / C +) (4-25) 

4+ = ±(1 — /”+)(2 + /”+) (4.26) 


Let 4? 5^,4, 4? an -d denote the expressions on the right sides of (1.35)— 
(1.40), in which 4 ^ + j 4^ - »"' are replaced by 4 ^ + j^n^~ 5 " From (1.35)-(1.37) and 
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(4.27) 


(1.38)— (1 -40), respectively, one obtains 

“S 1/2 = |N 1 ) +4 1 , +4 1) ] 
k + )$ 1/2 = | { [4° - (4i )+ «© i/a ] /( 2 »s )+ - i © i/2 

- [4 1 * - (4 ,+ «)Jt' /2 ] /(24 )+ - 1)m ,/2 } (4-28) 

« + )”t 1/2 = 5 { [4 i; - (4 )+ «© 1/2 ] /( 2 »S )+ - i)”t ‘ /2 

- [4 2> - (4 )+ «gt 1/2 ! /(24 ,+ - i)”t 1/2 } ( 4 - 2S >) 

for any ( j , k,n + 1/2) € f2i, and 

«#' = i [4 2) + 4 2) + 4 2) ] (4.30) 

(“DS 1 = 5 { [4 2) - (4 ,+ «)$'] /(24 )+ - 1 )S* 

- [4 2) - /(24 H - l© 1 } (4.31) 

K + )S’ = |{[4 2 ) -(4 h «)S 1 ]/(24 h -i)S 1 

- [4 2) - (4 H «)”,f ] /(2s-i? )+ - 1)$'} (4.32) 

for any ( j , &, n + 1) € fi 2 - 

Equations (2.18)-(2.56) are valid in the weighted-average v-t scheme. The accuracy of 
v-e scheme is also second-order in time and space. The stability condition is the same as 
that of a-e scheme with and v n replaced by / c+ and jP ?+ which are not constants any 
more. 

Test Problems and Discussions 

All of the two-dimensional schemes described in the previous section are tested us- 
ing standard model problems. The numerical results are compared with either the exact 
solutions or the numerical solutions obtained by some traditional finite-difference or finite- 
element methods to examine their accuracy, simplicity and effectiveness. 


15 



1. Test of the a-fi Scheme 

The translating Gaussian hill problem that was solved by using Taylor- Galerkin finite- 
element method in [5] is investigated here again. The governing equation is 


du du du d 2 u d 2 u . 

m + a ‘di + = "fe + W )y 


-1 < x < 1, -1 < y < 1 


which describes the problem that a Gaussian hill translates with a uniform velocity a = 
( a x ,a y ) and spreads isotropically with a diffusivity fi. Its exact solution is 


u(x,y,t) = -rr-rexp 


<r 2 (f) [ 2<r 2 (t) 


((* - *o - <M ) 2 + (y - yo - a y t?) 


where <x(t) = cr 0 y 1 + ifit/crfi, a x = acos<f>, a y = asin^>, in which <t> is the angle between 
the advection velocity and the x axis. In this problem, the initial conditions 


u(x,y, 0) 
u x (x,y, 0) 
Uy(x,y, 0 ) 


■ exp 


~2 (( x - x °) 2 + (y- yo?) 


2 cr§ 


tt(z,j/,0) 


^0 


(y - yo) 


are set as the exact solution at t=0, where cr 0 = 0.0707, <f> = 45°, x 0 = yo — —0.5 are used. 

The grid size used in the a-/x scheme is initially 25 x 86. Figure 3 shows the numerical 
solution and the distribution of the absolute error, which is defined as — u e |, for 
H — 0.0125 at t=0.5. The relative 1 l 2 error, which is defined as 


r _ VSgrV - ^)V(« - l)(m - 1) 

M = 

Umax 


obtained with At = 0.02 is 6.997 x 10 -4 and the numerical maximum value is 99.21% of the 
corresponding exact value. When the grid size is doubled to 50 x 173 and with At = 0.005, 
the relative Z -2 error is reduced to 3.84 x 10 -5 and the numerical maximum value becomes 
99.97% of the exact value. 

The same problem is then solved but with a much higher diffusivity of y, = 0.0625. 
Figure 4 shows the numerical solution and the absolute error distribution at t=0.5. The 
relative L 2 error of the numerical solution obtained on a 25 x 86 grid with At = 0.02 is 
1.788 x 10 -2 and the numerical maximum value is 96.07% of the exact value. Similiarly, 
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the relative L 2 error can be improved to 3.363 x 10 -3 and the numerical maximum value 
is up to 98.68% of the exact value by doubling the grid size to 50 x 173 with At = 0.005. 

Numerical results of the a-fj, scheme depict that its accuracy can be improved by using 
finer grid meshes and it can generate more accurate numerical solutions for lower diffusivi- 
ties. Once the grid size is fixed, the time step size will be the only parameter to determine 
the accuracy and stability. In the first test case with fi = 0.0125, different time step sizes 
have been used to obtain the best results as well as to satisfy the stability condition. In the 
second test case with v — 0.0625, the same time step sizes as those in the first case are used 
for a fair comparison between the numerical solutions for lower and higher diffusivity. The 
comparisons show that the numerical solutions for the lower diffusivity is more accurate. 


2. Test of the a-e Scheme 


The a-e scheme expressed in (2.12), (2.21), (2.22) and (2.15), (2.28), (2.29) is used to 
solve the problem that a Gaussian hill travels with a uniform advection velocity a = ( a x , a y ) 
without diffusivity, which is a special case of the previous Gaussian hill test problem. The 
governing equation is the pure advection equation 


du du du 
— + ax — + ay — 


0, -1 < x < 1, 


-1 < y < 1 


The exact solution and initial conditions are the same as those in the previous problem 
by letting /u, = 0. The exact solution is shown in Fig. 5(a). Numerical solutions at t=1.0 
obtained by using the same time step size At = 0.02 on a 50 x 173 grid are shown in 
Figs. 5(b)-5(d) for e = 0, 0.1, and 0.3, respectively. There are some oscillations along 
the direction of advection velocity when e = 0; they disappear when e = 0.1 or 0.3. But 
the relative L 2 error are 4.514 x 10 -3 , 2.83 x 10 -3 , and 3.08 x 10~ 3 , and the numerical 
maximum values are 98.2%, 97.6%, and 97.4% of the exact value for e — 0, 0.1, and 0.3, 
respectively. It can be seen that the numerical oscillation can be suppressed by increasing 
e and, meanwhile, the maximum value deviates more from the exact value because of the 
increased numerical diffusion. A large e may lead to too much numerical diffusion which is 
not expected. The solution with e = 0.1 is the best for this problem. 


3. Test of the // Scheme 


The n scheme is used to solve the unsteady diffusion equation 


du 

dt 


,d 2 u d 2 u . 


0 < y < 1 
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which was solved in [6, pp. 259-266] by using Crank-Nicolson implicit finite-difference and 
finite-element methods. 

The initial conditions are stated as follows: 


u(x , y, 0) = 20 + 80 [y — sin(0.57rx) sin(0.57ry)] 
u x (x , y , 0) = — 40?r cos(0.57tx) sin(0.57ry) 
u y (x,y,0) = 80 — 407r sin(0.57rx) cos(0.57ry) 

which are defined by the exact solution 

u(x, y , t) = 20 + 80 [y — e -0-5/i7r2 * sin(0.57rx) sin(0.57ry)] 

at t=0. The boundary conditions are: 


and 


u(0,y,t) = 20 + 80y 

u(l,y,t) = 20 + 80 

y — e °- 5fiw2t sin(0.57ry)j 

o 

<M 

II 

o' 

3 


u(x, l,t) = 20 + 80 

1 — e -0-5 ^ 2 * sin(0.5?rx)j 


%(0,y,t) = BO 

u y (l,y,t) = 80 [l — 0.57re -O-5/i,r2t cos(0.57rt/)| 

u x (x,0,t) = 0 

u x (x,l,t) = — 407re _O-5 ^ 7r2 * cos(0.57rx) 


Equations (3.4)-(3.9) are used at the boundary CEs. 

A 6 x 21 grid is used here for computation. Figure 6 shows the numerical solution at 
t=24000. The L 2 error is defined as 




T&r'jU'-uSj) 2 

\ (n — 1) x (m - 1) 


which is the same as that in [6]. The L 2 error is 5.9889 x 10~ 4 obtained by the y scheme with 
At = 100. Considering the results obtained by all of the numerical methods mentioned in 
[6] having the same accuracy (At, (Ax) 2 , (Ay) 2 ) as the y scheme, the best one is the Crank- 
Nicolson implicit finite-element method, whose errors are respectively L 2 = 6.86 x 10 -3 with 
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a 10 x 10 grid, and 1.59 x 10 -3 with a 20 x 20 grid. The numerical results can be improved 
so that those errors are reduced to 1.34 x 10 -4 and 8.2 x 10 -5 , respectively, by using the 
higher-order Crank-Nicolson method which has the ((At) 2 , (Ax) 4 , (Ay) 4 ) truncation error. 
By using the n scheme with more grid points, such as 30 x 60, an improvement is obtained 
having L 2 = 1.3467 x 10 -4 with At = 12. However, additional improvements cannot be 
achieved by further refinement of the grid, showing the need of a higher order accurate 
implicit scheme based on the CE/SE method. 


4. Test of the u-e Scheme 

Consider the inviscid Burgers equation 


l + S + | = 0 ' OS'* 1 ’ °*» S1 - 


where / = -it 2 , that has a discontinuous initial condition described by 


u(x,y,0) = 

1 

2 

if x < ^ 
“ 2 

and 

y < — 
y ~ 2 

= 

1 

1 

1 fx > - 

and 

1 

V> 2 

= 

1 

4 

otherwise 



for the first case, and 





u(x,y,0) = 

1 

2 

if x < ^ 

and 

y<l 
y - 2 

= 

1 

8 

lix> 2 

and 

1 

y> 2 

— 

1 

4 

otherwise 




for the second case. 

This problem contains complicated shock interaction, whose exact solution is referred 
to [7]. MacCormack’s scheme is used to solve the same problem to test the techniques of 
using a moving finer mesh which follows the shock as described in [7], where the fine grid 
point is up to 100 x 100, and the boundary conditions are specified as the exact solution 
for the first case, whereas for the second case, the conditions at the inflow boundaries at 
x = 0 and y — 0 are specified as the exact solution, but those at the outflow boundaries 
at x = 1 and y = 1 are computed by interpolation from the upstream grid points. For 
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the first case at t = 1.125, the exact solution is plotted in Fig. 7(a) and the numerical 
solutions obtained by the weighted-average v-t scheme with grid sizes of 50 x 173 and 
25 x 87 are shown in Figs. 7(b) and 7(c), respectively. Figures 8(a)-8(d) show the exact 
and numerical results for the second case at t = 1.125. In comparison with the numerical 
results by MacCormack’s scheme shown in Figs. 2 and 5 of [7], the present results in both 
cases are much better except a little smear near one of the discontinuities, yet the scheme 
is so simple and no additional ad hoc techniques are needed. Even with the coarse 25 x 87 
grid, the result is already much better. For a larger a , the numerical diffusion is larger and 
the smear near the discontinuity is more serious as depicted in Figs. 8(b) and 8(c). 

Conclusion 

The two-dimensional a-fj, and a-e schemes based on the space-time conservation element 
and solution element method have been used to solve the unsteady ad vection- diffusion equa- 
tion, unsteady pure advection equation, and unsteady diffusion equation. The a-e scheme 
has been modified to become the v-e scheme for solving the inviscid Burgers equation. 
Six numerical tests have been used to demonstrate that all of those schemes are superior 
to some of the traditional methods in accuracy and simplicity. The schemes themselves 
are physically sound, natural and simple. Highly accurate solutions can thus be obtained 
without using any ad hoc techniques. For further investigations, a more efficient implicit 
scheme with a higher accuracy based on the CE/SE method for solving the unsteady dif- 
fusion equation will be developed. 
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1/3 2/3 1 4/3 5/3 2 7/3 8/3 3 


° at time levels n=0, 1, 2, 3 ... 

• at time levels n=l/2, 3/2, 5/2 .... 


(a) Spatial positions of the mesh points marked 
by solid and harrow circles 


Fig. 1 The 2-D SEs and CEs in the CE/SE method 
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D' 



G’ = (j, k, n+1/2) 

A = (j+1/3, k+1/3, n) 

C = (j-2/3, k+1/3, n) 

E = (j+1/3, k-2/3,n) 

CE (1) (i, k, n+l/2)=box ABGFA'B'G’F' 
CE (1) (i, k, n+l/2)=box BCDGB'C’D'G' 

CE (I> (j, k, n+l/2)=box DEFGD'E'F'G' 

, k, n+1/2), i=l, 2, 3, and j, k, n=0, ±1, ±2, ... 


Fig. 1 Continued. 






G' = (j, k, n+1/2) 

B’ = (j-1/3, k+2/3, n+1/2) 

D' = (j-1/3, k-1/3, n+1/2) 

F = (j+2/3, k-1/3, n+1/2) 

SE 0) (j» k, n+1/2) = The union of four 
planes A’B’C'D'E’F, GBB"G", 
GDD"G" and GG"F'F and their 
immediate neighborhoods. 


(c) Solution element SE ( ) (j, k, n+1/2), j, k, n = 0,±1,±2,... 


Fig. 1 Continued. 



D' 



Fig. 1 Continued. 



E" 



Fig. 1 Continued. 



(b-w/3, h) 
/G~ 


(b+w/3, h) 
“B\ 


( 0 , 0 ) 

G 


>(2w/3, 0) 


(-b-w/3, -h) (-b+w/3, -h) 


0-2/3, k+1/3) (H /3 > k+2/3) ^ ^ 

"" :0.k) 


,.F 


0+1/3, k-2/3) 0+2/3, k- 1/3)' 


(f) The hexagon ABCDEF, the coordinate (£ XI ) and 
the parameters w, b and h. 


Fig. 1 Concluded. 
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j=1 j=2 j=3 j=nx 

nx is an even integer. ny is an odd integer. 

° at time levels n=2, 4, 6,.... 

• at time levels n=l , 3, 5,.... 

(a) The relative position of the mesh points 


0-2/3, 



0+1/3, 0+2/3, 

k-2/3) k-1/3) 


(b) Projected view of the interior CEs 

Fig. 2 The grid point distribution, the interior and boundary 
CEs used for solving the diffusion equation 
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(c) Projected view of the boundary CEs at x=0 and x=l 


(j-2/3, 

k+1/3) 




0 - 1 / 3 , 

k+2/3) 


(d) Projected view of the boundary CEs at y=0 



k-1/3) k-2/3) 


(e) Projected view of the boundary CEs at y=l 


Fig. 2 Concluded. 


29 




Fig. 3 The CE/SE solution and error distribution^ = 
50 ,m = 173, At = 0.005) at t=0.5 for v = 0.0125 
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Fig. 4 The CE/SE solution and error distribution^ = 
50, m = 173, At = 0.005) at t=0.5 for v = 0.0625 
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Fig. 5(a) The exact solution of purely advecting Gaus- 
sian hill at t=1.0. 
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Fig. 5(b) The CE/SE solution(n = 50, m = 173, e = 
0, A t = 0.02) of purely advecting Gaussian hill at t=1.0. 
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Fig. 5(c) The CE/SE solution(n = 50, m = 173, e = 
0.1, At = 0.02) of purely advecting Gaussian hill at 
t=1.0. 
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Fig. 5(d) The CE/SE solution(n = 50, m = 173, e = 
0.3, At = 0.02) of purely advecting Gaussian hill at 



Fig. 6 The CE/SE solution(n = 6, m = 21, At = 100) of 
diffusion equation at t=24000. 


36 




Fig. 7(a) The exact solution of the inviscid Burgers 
equation at t= 1.125 for the first case. 
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Fig. 7(c) The CE/SE solution (n = 25, m = 87, e = 
0.5, a = 1, At = 0.01) of the inviscid Burgers equation 
at t= 1.125 for the first case. 
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Fig. 8(a) The exact solution of the inviscid Burgers 
equation at t=1.125 for the second case. 
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Fig. 8(c) The CE/SE solution (n = 50, m — 173, e = 
0.5, a — 1 , At = 0.01) of the inviscid Burgers equation 
at t=1.125 for the second case. 
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Fig. 8(d) The CE/SE solution (n = 25, m = 87, e = 
0.5, a = 1, At = 0.01) of the inviscid Burgers equation 
at t=1.125 for the second case. 
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